A lattice mesoscopic model of dynamically heterogeneous fluids 
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We introduce a mesoscopic three-dimensional Lattice Boltzmann Model which attempts to mimick 
the physical features associated with cage eflects in dynamically heterogeneous fluids. To this 
purpose, we extend the standard Lattice Boltzmann dynamics with self-consistent constraints based 
on the non-local density of the surrounding fluid. The resulting dynamics exhibits typical features of 
dynamic heterogeneous fluids, such as non-Gaussian density distributions and long-time relaxation. 
Due to its intrinsically parallel dynamics, and absence of statistical noise, the method is expected 
to compute significantly faster than molecular dynamics, Monte Carlo and lattice glass models. 
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Very slow relaxation and long equilibration times are 
a typical feature of complex fluids, such as supercooled 
liquids, polymers and others [EQl- The numerical study 
of these complex systems is usually undertaken by many- 
body simulation methods, such as Molecular Dynamics 

tjMonte Carlo 0, and various types of lattice 'glasses' 
Ig, Q) B S ll^ ■ Since many-body simulations are com- 
putationally demanding, it is worth exploring whether 
the salient features of dynamic heterogeneities can be 
reasonably described by effective one-body techniques in 
the spirit of density functional theory A power- 

ful one-body technique is the Lattice Boltzmann (LB) 
method. The LB dynamics consists of three basic steps: 
Free-streaming, collisional relaxation ^.nd (effective) 
intermolecular interactions |l3 . Judicious choice of these 
intermolecular interactions permits to describe the dy- 
namics of a variety of complex flows However, appli- 
cability of LB to glassy-like fluids remains an open prob- 
lem |l5| . A crucial aspect of the physics of complex fluids 
are dynamic heterogeneities and geometrical frustration, 
i.e. the presence of sterical constraints which reduce the 
phase-space available to the fluid system. Real systems 
like colloids, or granular materials, exhibit glassy dynam- 
ics associated to jamming: Density, rather than temper- 
ature, is the dominant effect in slowing down strongly 
the dynamics In an attempt to model these ef- 

fects, we develop a lattice Boltzmann equation (LBE) 
in which free-particle motion is confined to a subset of 
links which fulfill self-consistent dynamic constraints on 
the surrounding fluid density. The presence of these ki- 
netic constraints on free particle motion leads to drastic 
departures from simple fluid behavior, such as long-time 
relaxation and non-Gaussian density fluctuations. We 
begin by considering a standard LBE with a single relax- 
ation time [l7llli|: 

Mr, t) - f,{r,,t - At) = -iuAt [/, - /f ] (r„ t - At) (1) 

where Ax and At are space- and time-steps, respectively. 



fi{r,t) = /(r,v = Ci,t) is a discrete distribution func- 
tion of particles at site r and at time t moving along 
the direction i of a lattice, with discrete speed c^, and 
r; = r — CiAt. The previous equation can be seen as 
the combination of collision and streaming steps. In the 
collision, the distribution functions fi relax to a local 
equilibrium ff in a time lapse of the order of uj^^, such 
that the distribution function after a collision /f is: 

f^{r,,t - At) ^Mr,,t~ At) - ujAt [/, - /f] {r,,t - At) 

(2) 

The distribution function /f then streams freely to neigh- 
bor sites, so that the updated values are just the shifted 
post-collisional distributions: 



/,(r,i) = /f(r„<-Ai) 



(3) 



The large-scale behavior of system depends crucially on 
the form of the local equilibria. In the present work, the 
equilibrium distribution functions ff are expressed as: 



/f(r,t) = Wip{r,t) 



(4) 



where Wi is a set of lattice-dependent weights normalized 
to unity. The local density p(r, t) in eq. Q is obtained 
by a direct summation upon all discrete distributions: 



p(r,t) =^/,(r,i). 



(5) 



Since the local equilibria do not depend on the local fluid 
speed, the only conserved quantity is the fluid density, 
which means that in the continuum limit, the system 
obeys a simple diffusion equation: 9tp(r, t) — DW^p{r, t) 

with diffusion constant D — (zJSt ~ Should 

flow phenomena be of interest as well, we should just add 
a quadratic term in the flow speed to the local equilibrium 
distribution functions In the following, we shall refer 
to a three-dimensional cubic lattice of size Lx Lx L with 
Co = (0,0,0) and a = (±Aa;/Ai, 0, 0),(0, ±Aa;/Ai, 0), 
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(0,0,±Aa;/At), i = 1,2,..., 6. In this case wq = 1/3 
and Wi = 1/9 for a = Ax/ At, i = 1,2,..., 6. We en- 
force sterical constraints which have proven to be ef- 
fective to capture the physics of glassy systems where 
density is the dominant observable [3, 13 • Kinetic con- 
straints on the evolution of the system are enforced by 
the following functional rule: Propagation from a site r 
to one of its 6 neighbors r' is permitted only if the non- 
local densities p„z(r) = /^("^ + c^At) and Pniiy') = 
Si=i p(r' + CiAt) prior and after streaming, respectively, 
both lie below a given density threshold, S. In the limit 
5 — i- cw, the effective propagator taking the system from 
site r at time t to site r' at time t+At, reduces to the stan- 
dard free-particle form, G'i(r,r'; At) = 5(r' — r — c^Ai). 
In the opposite limit, S" — > 0, no motion is allowed and 
Gi — > b(y' — v) at all sites, corresponding to structural ar- 
rest. It is therefore clear that the ratio < p > / S {< ... > 
stands for a spatial average), serves as a control param- 
eter driving the system from the purely diffusive to the 
structural arrest regime. The above rule is implemented 
as follows: 1. Initialize the system by randomly choosing 
N lattice sites and set them at a local density po. The 
remaining — N sites are set at density 0. The aver- 
age density in the system is then < p >= xPo < Pm, 
where pM = Po is the maximum possible average den- 
sity in the system, and x = N/ is the concentration of 
'loaded' sites; 2. Compute local densities p(r) via eq. ijSl; 
3. Compute the equilibrium distribution functions /f via 
eq. ^. Perform the collision ||2J) on all the lattice sites 
to compute ; 5. Look for all the lattice sites r* such 
that P{^* + Ci) < S where S" is a fixed threshold; 
6. Perform a pre-streaming, according to eq. of the 
ff computed at step 4 only along links emanating from 
the lattice sites r* and pointing towards neighbor sites 
r* ; 7. Compute again local densities p* (r) on all the lat- 
tices sites; 8. Look for all the lattice sites r** such that 
Xli=i + Ci) < S] 9. Perform the effective stream- 
ing step of the ff computed at step 4 only along links 
from r* to r** sites (these links will be denoted as active 
ones), otherwise the ff are not moved; 10. Go to 2.. The 
LB scheme described above is intrinsically distinct from 
those used for non-ideal fluids. Indeed, while the latter 
include potential energy via effective interactions which 
leave the free-propagator (kinetic energy) unaffected, in 
our case the interactions alter the structure of the kinetic 
energy operator. 

We have simulated two lattice sizes, L = 16 and 32, 
with Aa; = At — 1. We changed lu in the range [0.1 : 1] 
and did not find any dependence of results on its spe- 
cific value. Thus, we set uj = 0.1. We used po = 0.5 
and S — 1.5, corresponding to single site density thresh- 
old ps — S/& = 0.25. Keeping po fixed, we chose the 
smallest value of S such as to ensure sluggish dynam- 
ics at high densities, while still allowing the system to 
be uniform at low densities. By running several simu- 



lations with Po = 0.5, we have found that for S* > 1 
the system evolves by diffusive smoothing of the density 
gradients towards a long-time state characterized by a 
uniform density pattern when < p 0.03. Moreover, 
when S ^ 2.3, the system does not show any singular 
behavior for densities smaller than the maximum pos- 
sible one, Pm (see below). This latter feature is also 
observed by increasing po while keeping S = 1.5. This 
is because, by increasing po at a fixed average density 

< p >, the number A'' of lattice sites to be initialized 
with density po decreases, and consequently it becomes 
more difficult for the kinetic constraints to be effective. 
In conclusion, the system seems to exhibit a non-smooth 
transition from diffusive to sluggish behavior as the re- 
duced density A =< p > / ps is increased. 

This is shown in Fig. 1, where we plot the order pa- 
rameter m = {pmax — Pmin)/ Po as a function of A by 
keeping < p >= 0.12 and po — 0.5 fixed and varying S. 
Here, Pmax and pmin are the maximum and minimum 
values of p, respectively, at steady-state. The values of 
TO were averaged over 50 independent runs on systems of 
size L = 32. The simulations yield to ~ for A ^ 0.36, 
indicating a purely diffusive behavior, then to undergoes 
a sharp rise and the system enters the sluggish regime. 
At values of A ^ 1.44, the system is nearly frozen in its 
initial configuration and the order parameter to settles 
down to 1 since pmax = Po and Pmin = 0. Figure 1 
clearly indicates the existence of three distinct regimes, 
namely a low-density diffusive regime at A < Ad — 0.36, 
a high-density frozen regime at A > Xp — 1.44, and a 
sluggish regime at intermediate densities Ajj < A < Af- 
Since the initial density of occupied sites is po, the initial 
non-local density Pniiy) = Y^^i=iP{^ + c^At) can be at 
most 6po. Therefore, the condition for a purely diffusive 
behavior is 6po ^ A < 0.24, which is in a reasonable 
agreement with the value provided by simulations. The 
value of A marking the transition from the glassy to the 
frozen regime can be determined by considering that the 
average initial non-local density is Pn/(r) = 6xPo- When 
6xpo ^ S => A > 1.0, the system is likely to be frozen. 
This underestimates the value obtained from simulations, 
but the value of the order parameter to(1) ~ 1.3, well be- 
low its maximum to ~ 2.4, indicates that the system is 
approaching the frozen regime. We choose A = 0.48, in 
order to select a regime with clear departure from ideal 
fluid behavior (to(0.48) ~ 2.35). For this set of param- 
eters the system evolves from the initial random con- 
figuration forming some clusters until, at long times, it 
gets arrested in one of the highly heterogeneous states. 
The arrest time decreases with increasing average density 

< p > . To analyze the qualitative difference with respect 
to the simple diffusion dynamics, we also inspected the 
probability distribution functions of density. It is ob- 
served that the system, which starts from a two-peak 
distribution at p = and po, rapidly fills- up all available 
values in the range [0 : 1.1]. At long times there is a sharp 
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FIG. 1: The order parameter m (•) as a function of the re- 
duced density A with po — 0.5, < p >= 0.12, and L = 32. 
Also shown is the steady-state value of participation number p 

(□), defined as p = 1/{L' Ef=iP?)' with = p(r,)/(xi'po) 
being the occupation probabihty of the site i. By definition 
p{t = 0) = X s-iid p = 1 in the case of uniform distribution. 
This plot shows that x < P < 1 indicating that the system 
never becomes more localized than in the initial configuration. 



peak at p ~ 0.05, which corresponds to the background 
density and the distribution is non uniform everywhere 
else, with a pronounced shoulder at p ~ 0.1. 
We computed the time autocorrelation function 



FIG. 2: Autocorrelation function h{t) as a function of time 
t for < p 0.08(A), 0.09(o), 0.10(*), 0.11(«), 0.12(*), 
0.18(0), 0.30(n) with po = 0.5, S = 1.5, and L = 32. 



h(t) 



« 5p{t + to)5p{to) » 
« 5p{tQ)dp{to) » 



(6) 



where << ... >> denotes an average over space and ini- 
tial times ^0 and Sp{t) — p{t)— < p >. The plots for 
several values of the initial density < p > are shown in 
Fig. 2 for the case L = 32. Data were obtained by aver- 
aging over 50 independent runs for each value of < p >. 
For very small values of < p >, the system goes to a final 
state with uniform density and h(t) relaxes to zero. By 
increasing < p >, h{t) starts forming a plateau and stays 
close to the unit value for a time span which increases 
rapidly with increasing mean density < p >. Even at 
high densities, the correlator does not show the "two- 
step" relaxation behavior often found in glassy materi- 
als. We believe that this is due the absence of a rattling 
motion in our model (a similar behavior is found in the 
Kob- Andersen (KA) model 0| for lattice glasses). 

We compared our results with the predictions of mode- 
couphng theory (MCT) 0. MCT predicts a power-law 
time decay away from the plateau. We therefore tried 
to fit the short-time behavior of the function h{t) with 
a power law of the form / — Bt'', where /, B and b 
are fitting parameters. Such a power-law decay is in- 
deed reproduced by our simulations in the high-density 
regime where from all the fits we find / very close to 1 
(to within 10~^) as in the KA model |5j. The coefficient 
B decreases at increasing values of < p > varying in the 



range [4 x 10"^ : 3 x 10""]. At variance with MCT pre- 
dictions, the power-law exponent b is not independent 
of density but decreases at increasing values of < p > 
varying in the range [0.8 : 1.0]. This density dependence 
of b is found also in the KA model Q with b decreas- 
ing at increasing values of < p > varying in the range 
[0.8 : 1.1], which is consistent with our results. In MCT 
the decay at longer times is predicted to be a stretched 
exponential, h(t) cx exp {—{t/r)^), where the exponent (3 
is density-independent and the relaxation time t is the 
relevant physical parameter. Indeed, this time scale in- 
creases strongly as density is increased (MCT predicts 
a singular behavior for t at a density smaller than the 
maximum one pm)- Chemically different materials re- 
lax in a qualitatively similar manner with relaxation 
functions obeying the Kohlraush- William- Watts function 
exp(— (i/r)^) j2flj. We fitted successfully h{t) at long 
times for < p >> 0.12 by using a stretched exponential 
with T and [3 as fitting parameters and found that the 
inverse relaxation time 1 /t vanishes at a critical density 
Pc = 0.412 ± 0.010, with power law 6(pc- < p >)^=, be- 
ing 7c = 4.68 (in the KA model the critical exponent 7c 
is ^ 5 0), while the stretching exponent (3 is density- 
dependent. MCT predicts such a power law behavior for 
the relaxation time with a system-dependent exponent 
7c 01 • It is remarkable that, despite the density depen- 
dence, we find numerically that /?(< p >) ~ 6(< p >), 
as predicted in MCT for the time dependence of the au- 
tocorrelation function j^J]. We stress that the critical 
value pc is smaller than the maximum density pM, cor- 
responding to a fully-loaded lattice. The plot of the rele- 
vant physical quantity 1/r as a function of pc— < p > is 
shown in Fig. 3, where we also report the fitting values 
of 1/r for the system size L = 16. Also in this case, we 
found that l/r vanishes with a power-law behavior and 
the estimated critical density is 0.405 ± 0.029, which is 
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consistent with pc within the error range, with no signifi- 
cant lattice size effect. This suggests a singular behavior 
at about < p >= 0.412. We also inspected the ratio 
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FIG. 3: Inverse relaxation time 1/t as a function of 0.412— 
< p > for po = 0.5, 5 = 1.5 and L = 16(A), 32(o). The 
full line has slope 4.68. In the inset the fraction F of active 
links versus 0.471— < p > is shown for the same parameters. 
The full line has slope 4.19. For a comparison we plotted the 
results of the KA model by using dashed lines. 



F of active lattice links to the total number 6L^ of lat- 
tice links. This quantity can be viewed as a measure of 
the degree of glassiness, as observed in kinetically con- 
strained lattice glass models since we expect the 
relaxation time r to diverge in the limit F —^ 0. In the 
limit <p>— i'Oji^^l since the kinetic constraints do 
not play any role. When < p >—> pc, we expect F —> 
since in this limit the kinetic constraints are very effective 
in slowing down the dynamics j^. 

In the inset of Fig. 3 we plot the steady-state values 
of F, averaged over 50 runs, for each value of < p >. 
We found that F vanishes at a critical density p^ — 

0.471 ± 0.005, with power law 120(p^- < p >)^'-, be- 
ing 7^ = 4.19. In the case of the KA model, a similar 
power law is found with critical exponent 7^ ^ 4.7 J^. 
In the same inset we also report the values of F for the 
system size L — 16. Also in this case we found that F 
vanishes with a power law behavior and the estimated 
critical density is 0.482 ± 0.016, which is consistent with 
p^ within the error range, and no significant lattice size 
effect. It is interesting to note that p^ ~ pc, supporting 
the conclusion that our model shows singular behavior at 
a density smaller than the maximum possible one, pm. 

Summarizing, we have introduced a mesoscopic LB 
model which appears to reproduce some physical features 
of dynamically heterogeneous fluids, such as sluggish re- 
laxation and continuum density distributions. To this 
purpose, the standard LB dynamics has been augmented 
with self-consistent constraints based on the non-local 



density of the surrounding fluid. A typical run on a 32"^ 
lattice (four times larger than typical lattice glass simu- 
lations [aLDl) takes just a few minutes on a 2.4 GHz Intel 
Xeon processor. 
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